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Abstract 

The  multigrid  waveform  relaxation  (WR)  algorithm  has  been  fairly  stud¬ 
ied  and  implemented  for  parabolic  equations.  It  has  been  found  that  the 
performance  of  the  multigrid  WR  method  for  a  parabolic  equation  is  practi¬ 
cally  the  same  as  that  of  multigrid  iteration  for  the  associated  steady  state 
elliptic  equation.  However,  the  properties  of  the  multigrid  WR  method  for 
hyperbolic  problems  are  relatively  unknown.  This  paper  studies  the  multi¬ 
grid  acceleration  to  the  WR  iteration  for  hyperbolic  problems,  with  a  focus 
on  the  convergence  comparison  between  the  multigrid  WR  iteration  and  the 
multigrid  iteration  for  the  corresponding  steady  state  equations.  Using  a 
Fourier-Laplace  analysis  in  two  case  studies,  it  is  found  that  the  multigrid 
performance  on  hyperbolic  problems  no  longer  shares  the  close  resemblance 
in  convergence  factors  between  the  WR  iteration  for  parabolic  equations  and 
the  iteration  for  the  associated  steady  state  equations. 
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1  Introduction 


The  advent  of  a  new  generation  of  massively  parallel  computers,  consisting 
of  hundreds  or  thousands  of  processors,  has  caused  previously  unattrac¬ 
tive  numerical  algorithms  to  be  reexamined.  For  a  numerical  algorithm 
to  fully  exploit  the  power  of  such  machines,  it  must  be  decomposable  into 
largely  independent  pieces  which  can  be  distributed  among  the  available 
computer  processors.  The  waveform  relaxation  (WR)  method  (or  dynamic 
iteration  method),  originally  proposed  for  VLSI-simulation  [8]  [12],  is  such 
a  method  for  solving  systems  of  ordinary  differential  equations  (ODEs)  and 
time-dependent  partial  differential  equations  (PDEs).  It  decomposes  a  full 
system  into  smaller  subsystems  which  can  be  solved  concurrently.  Addition¬ 
ally,  it  allows  different  integration  step  sizes  to  be  used  for  different  subsys¬ 
tems,  resulting  in  substantial  savings  in  computation  for  some  applications 
[7]  [10]  [13]  [14]  [15]  [19]. 

Many  studies  on  the  WR  method  and  its  acceleration  techniques  have 
been  made  on  parabolic  problems  [4]  [6]  [9]  [16]  [17].  In  those  studies,  an  im¬ 
portant  approach  has  been  to  establish  a  quantitative  comparison  between 
the  WR  iteration  for  time-dependent  PDEs  and  the  associated  iteration  for 
the  corresponding  steady  state  problems,  called  static  iteration,  since  the 
latter  has  been  extensively  investigated.  The  studies  have  found  that,  for 
parabolic  problems,  the  convergence  rates  of  the  WR  iteration  and  static 
iteration  are  quantitatively  quite  close  [11].  When  a  multigrid  technique  (in 
space)  is  incorporated  into  the  WR  method,  the  typical  multigrid  accelera¬ 
tion  can  be  achieved  with  a  rate  that  is  a  small  perturbation  from  the  one 
for  the  associated  static  multigrid  iteration  [9]  [16].  Therefore,  the  result¬ 
ing  multigrid  WR  method  significantly  increases  the  speed  of  convergence 
in  solving  parabolic  problems,  making  it  competitive  with  the  traditional 
time-stepping  methods. 

For  hyperbolic  systems,  the  WR  method  has  been  shown  to  be  concep¬ 
tually  similar  to  the  semi-discrete  subdomain  iteration  and  share  the  same 
convergence  properties  as  for  parabolic  problems  [1].  Therefore,  it  is  natural 
to  consider  a  multigrid  acceleration  to  the  WR  iteration  because  the  tech¬ 
nique  has  been  effectively  used  for  solving  non-eUiptic  steady  state  problems 
[3].  For  certain  classes  of  non-eUiptic  steady  state  problems,  the  multigrid 
technique  has  substantially  improved  the  speed  of  convergence.  Since  it 
can  be  adapted  into  the  WR  method  for  solving  time-dependent  hyperbolic 
PDEs  in  the  exact  same  format,  a  question  is  raised  about  whether  its  per¬ 
formance  is  stiU  analogous  to  its  steady  state  counterpart  as  reported  by 


1 


Brandt  [3]. 

This  paper  addresses  the  issue  of  performance  differences  of  the  multigrid 
WR  method  on  parabolic  and  hyperbolic  problems.  It  studies  the  conver¬ 
gence  factors  of  the  multigrid  WR  method  for  problems  whose  space  operator 
is  not  elliptic,  or  in  which  the  elliptic  principal  part  of  its  space  operator  is 
small.  Many  problems  in  fluid  dynamics,  and  in  other  fields,  are  of  these 
types.  The  Fourier  method  was  shown  to  be  a  powerful  tool  for  convergence 
analysis  of  numerical  schemes  for  these  types  of  problems  [3].  In  this  paper, 
a  Fourier-Laplace  analysis,  i.e.,  Fourier  analysis  in  space  and  Laplace  trans¬ 
form  in  time,  is  used  in  two  case  studies  that  involve  hyperbolic  equations. 
Interesting  phenomena  were  observed.  In  particular,  it  was  found  that  for  a 
hyperbolic  equation,  the  multigrid  WR  iteration  could  fail  to  converge  while 
its  static  analogue  works  well  on  the  associated  steady  state  problem. 


2  Fourier-Laplace  Analysis 

The  Fourier  method  is  an  indispensable  tool  for  analyzing  both  differential 
equations  and  discrete  solution  methods  for  time-dependent  problems.  It  is 
especially  powerful  for  the  convergence  analysis  of  multigrid  iteration,  be¬ 
cause  it  provides  insight  into  the  details  of  the  basic  interaction  between 
the  coarse  grid  correction  and  fine  grid  relaxations.  For  steady  state  model 
problems,  comprehensive  Fourier  analysis  has  been  developed  [2]  [18].  For  a 
parabolic  model  problem,  the  practical  convergence  estimates  of  the  multi¬ 
grid  WR  iteration  were  obtained  by  a  Fourier-Laplace  analysis  [16].  In  this 
paper,  such  analysis  is  used  in  two  case  studies  that  involve  hyperbolic  equa- 
tions  formulated  as 

^  +  Lu  =  f,  t>0,  =  0)  =  uo,  (1) 

at 

with  periodic  boundary  conditions.  L  is  an  m-dimensional  linear  operator 
with  constant  coefficients  defined  in  an  infinite  space.  In  order  to  employ  a 
multigrid  WR  iteration  to  (1),  the  equation  is  first  discretized  in  space 

^  +  LHUh  =  fh,  t>0,  Uh{0)  =  uo,  (2) 

at 

where  Lh,  Uh,  and  fh  are  discrete  approximations  to  X,  u,  and  /  obtained 
by  spatial  finite  differences.  Then  the  multigrid  WR  iteration  is  applied  to 
(2)  with  h  as  the  finest  grid  size.  Our  analysis  is  restricted  to  a  two  grid 
WR  cycle  described  as  follows. 
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1.  Let  Lh  be  split  a.s  Lh  =  M  -  N.  Perform  vi  pre-smoothing  steps: 


=  Nui;~^'>+fh,  t>0,  ul"^(0)  =  «o,  v  =  l,2,---,vi 

where  the  starting  function  is  given. 

2.  Restrict  the  defect  from  grid  h  to  grid  H : 

J  (^l) 

+  dH-=li^dh.  (3) 

3.  On  the  coarse  grid,  solve 

-f-  Lhwh  =  dn,  Wi?(0)  =  0. 

4.  Correct 

Uh  -  -  Ih'^H 

where  is  a  suitable  interpolation  from  grid  H  to  grid  h. 

5.  Perform  V2  post-smoothing  steps  on  Uh- 

The  error  of  a  complete  two  grid  WR  cycle  described  above  satisfies 

4'^  =  4\^)  =  0,  (4) 

where  V  is  the  two  grid  WR  iteration  operator.  The  Laplace  transform  of 


(4)  is 

(5) 

with 

V{z)  =  S''^{z)CG{z)S’^^{z), 

(6) 

and 

CG{z)  =  1-  I^{z  -h  Lh)-^I^{z  -I-  Lh). 

(7) 

The  matrix  functions  CG{-)  and  5(-)  are  the  Laplace  transforms  of  convo¬ 
lution  kernels  of  the  coarse  grid  corrector  and  the  smoother  respectively. 
Note  that  CG{0),  5(0)  and  y(0)  are  respectively,  the  coarse  grid  correction, 
smoothing  and  two  grid  iteration  operators  for  the  corresponding  steady 
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state  problem  LhUh  =  A*  A  detailed  derivation  of  (6)  and  (7)  can  be  found 
in  [9]. 

Assuming  that  all  the  entries  of  are  rational  functions  of  z  vanishing 
at  infinity  with  poles  having  negative  real  parts,  and  taking  V  as  an  operator 
on  i/P(i2+,C”)  (1  <  p  <  oo),  the  spectral  radius  p(V)  =  limfc^oo 
satisfies  (see  [9]) 

p{V)  =  max  p{V{z))  =  mdi.x  p{V{z)).  (8) 

Since  the  multigrid  technique  is  used  in  space  only,  the  convergence  estimate 
of  the  multigrid  WR  iteration  can  be  obtained  by  performing  a  Fourier 
analysis  on  each  of  equations  (5). 

Let  1  •  1  denote  the  max-norm  in  (note,  m  is  the  dimension  of  the 
PDE  (1)).  Define  the  frequency 

where 

|^^|<|,  l^^'l  <7rand0j’ +  i  =  2,---,2’", 

in  which,  W  :=  iriiuh,  •  •  • ,  «m),  4  =  0  or  1,  and  at  least  one  of  4’s  equals 
1.  The  exponential  Fourier  mode  on  h  grid 

exp{i0^ xlh)i  X  €  C™, 

is  an  infinite  dimensional  vector  determined  by  the  grid  points,  and  appears 
as  the  mode  €xp{i20h/H)  on  H  =  2h  grid.  Therefore,  on  H  grid,  it  co¬ 
incides  with  aU  Fourier  modes  exp{iO^  x/h),  j'  =  1, •  •  •,2’”,  j'  7^  A  Thus, 
restriction  operators  introduce  coupling  between  each  lower  mode  9^  and  its 
(2”^  - 1)  high  frequency  harmonics  9^ ,  j  =  2,- ■ -,2^.  Interpolation  operators 
introduce  coupling  among  the  same  modes.  This  can  be  represented  as 

I^exp{i9^xlh)  =  ih{9^)exp(i29^xlH),  y  =  l,-*-,2™,  (9) 

and 

I^€xp{i29h/H)  =  ^  i^{9^)exp{i9^x/h).  (10) 

If  the  set  of  Fourier  modes  in  h  and  H  grid  are  denoted  by 

Xh{9)  =  [exp{ie^xlh),  ■  ■■,exp{i9‘^"xlh)]  and  Xh{9'^)  =  [exp{i29^xl H)], 
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the  equation  (9)-(10)  can  be  written  as 

i^Xh(0)  =  i^{0)  = 


and 


iifX„ie^)  =  Xh{0)iUff),  1^(0)  = 


The  matrices  I^i^)  and  1^(0)  are  called  the  matrix  symbols  of  and 
respectively.  In  this  paper,  the  symbol  of  a  matrix  A  will  be  denoted  by 
A{0). 

If  the  operators  Lh  a.nd  Lh,  and  the  smoother  in  Laplace  domain  S{z)  do 
not  introduce  coupling  of  more  Fourier  modes,  then  the  set  of  2"*  harmonic 
modes  Xk{0)  is  an  invariant  subspace  of  iteration  operator  V{z)  satisfying 


Viz)Xh{0)  =  Xh{0)Vi9,  z)  for  aU  0, 


where  the  matrix  symbol  of  V(z)  has  the  form 

V(0,z)  =  Jiez  >  0. 

The  spectral  radius  of  F(z)  can  then  be  obtained  by  collecting  its  value  on 
each  set  of  2^  Fourier  modes 


p(V(z))  =  snpp(V(0,z))=  max  p(V(0,z)). 

0  1^  l<f 

Combining  (8),  the  asymptotic  convergence  rate  of  the  two  grid  WR  iteration 
can  be  calculated  by 


p(V)  =  max  max  p(V(6,z)) 
^  Rez=0\d^<Y^  ^ 


(11) 


3  Case  Studies 

In  order  to  assess  the  multigrid  performance  on  hyperbolic  equations,  it  is 
useful  to  start  from  the  case  study  on  the  model  hyperbolic  equation 


Ui  +  aUx  +  bUy  =  /.  (12) 

A  wide  class  of  discretization  methods  for  equations  of  the  above  form,  even 
in  the  context  of  a  system  of  PDFs,  involves  central  differencing  of  the 
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original  equation  with  an  additional  term  which  is  used  for  stabilization. 
This  additional  term,  called  artificial  viscosity,  may  be  of  different  order 
and  usuaUy  corresponds  to  a  discretization  of  elliptic  operator  times  a  smaU 
coefficient.  Here,  two  such  discretization  schemes  are  chosen  for  the  case 
studies.  The  analysis  of  these  schemes  are  done  using  First  Differential 
Approximation  (FDA)  [-3]. 


Case  1.  Consider  the  equation 

ut  +  aux  +  huy  -  =  / 

with  given  initial  condition  and  periodic  boundary  condition.  This  is  the 
FDA  for  the  following  discretization  of  (12): 

1  r  V2-73 

Lh-.=  T  4/3  a/2-/3  . 

^  [  -6/2  -  /? 

Note,  Lh  is  a  first  order  discretization  to  a-§^  +  6^.  The  matrix  symbol  of 


Lh{0^) 


Lh{9)  = 


,  e=  =  {9{,e{),  j  =  l,---,4,  (13) 


Lh{9^) 


A9{si'n?{e{l2)  +  sin2(^/2))  +  i{asind{  +  bsin9l) 

Lh{0^)  =  —  ^  • 

The  matrix  symbol  of  the  coarse  grid  operator  is 

i/f(2^)  =  [M20')]-  (14) 

Let  the  bi-linear  interpolation  be  chosen  for  and  the  restriction  operator 
be  =  {IhY-  Their  matrix  symbols  are  given  by 


y  =  l,--.,4. 
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Symbols  for  the  damped  Jacobi  and  red-black  Gauss-Seidel  WR  smoother 
are  found  to  be 


Sj{0,z) 


sie\z) 


s{e\z)  \ 


(15) 


S{0\z)=^-^{d-u;Lh{0n),  0<a;<l,  i  = 

ZLO  “r  w 


Sa{0,z) 


A 


Sgs{0,z)  = 

1  +  A  “(1  A  A) 
1-A  -(1-A) 


SAiO,z)  0 
0  SB{e,z) 

,  SBi0,z)  =  ^ 


1  +  B  -{1  +  B) 
1-B  -(1-5) 


A: 

B  = 


_  2^{cos0l  +  cosBl)  -  i{asinO\  +  bsin6\) 
4/3  “}-  zh 

2l5{—cos9\  +  cosOl)  —  i{—asin9l  +  hsin9\) 
40  +  zh 


Tables  1  and  2  list  the  computed  spectral  radii  of  the  two  grid  WR 
operators  for  the  damped  Jacobi  relaxation  and  Red-Black  Gauss-Seidel 
relaxation  based  on  (11).  The  number  of  smoothing  steps  is  chosen  as  u  = 
vi-\-V2  =  1.  For  a  comparison,  using  the  same  number  of  sample  points  for 
the  space  variables,  the  computed  spectral  radii  of  the  related  two  grid  static 
iteration  operators  are  listed  inside  parentheses.  These  data  show  that  the 
performance  of  the  two  grid  WR  iteration  on  this  problem  is  quantitatively 
quite  close  to  its  static  analogue,  a  phenomenon  that  has  been  observed  and 
theoretically  justified  for  parabolic  problems. 


Table  1:  Spectral  Radius  of  Two  Grid  Damped  Jacobi  WR  (w  =  2/3,  6=1) 


a 

= 

^=.5 

^  =  .75 

/?  =  1. 

1.0 

.8416  (.8357) 

.7086  (.7001) 

.6072  (.5816) 

.5887  (.5644) 

0.5 

.7644  (.7480) 

.7509  (.7151) 

.7815  (.7298) 

.7992  (.7789) 

0.0 

1.0  (1.0) 

.9999  (.9999) 

.9999  (.9999) 

.9999  (.9999) 
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Case  2.  Next,  consider  a  third  order  approximation  of  (12),  for  which 
the  FDA  is 

Ut  +  au^  +  buy  +  h^PA^u  =  /, 

i.e.,  the  operator  +  6^  is  approximated  by 


Lh 


P 

2P  -8P  +  6/2  2p 

P  -8^ -a/2  20/?  -8/?  +  a/2  P 

2p  -SP  -  6/2  2P 

P 


The  matrix  symbols  for  Lh  and  Lh  are  again  represented  by  (13)  and  (14), 
but  with  the  elements 


.  16P(sirP(ei/2)  +  sin‘^{ey2)f  +  ijasinej  +  bsinOj) 

Lh{^  )  —  ~  ""  T”  ’ 


The  damped  Jacobi  WR  smoother  has  the  matrix  symbol  (15)  with 

S{e^,z)  =  -^id-u:Lh{e%  d=^,  0<a;<l,  J  =  l,---,4. 

Thus,  it  has  exactly  the  same  Fourier  smoothing  factor  (see  [18])  as  its  static 
analogue 

pr  =  max  max  |5(0-',2:)|  =  max  |5(0-’,O)| 

=  max  |1  -  oj{-{sinP{d{j2)  +  sin^{9{l2))^  +  r— (asira^  +  6sm^))|. 
2<j<4  5  2U/? 

Let  w  =  1/2.  Since  7r/2  <  <  tt,  j  =  2,3,4,  the  Fourier  Jacobi  WR 

smoothing  factor  satisfies 


pT  ^ 


|a|  +  1^1 

40/? 


)2<  1 


when  P  > 


(p  +  1^1 

32 


This  implies  that  the  Jacobi  WR  smoother  is  as  effective  as  the  associated 
static  smoother  in  eliminating  the  error  with  oscillatory  modes. 

Using  the  same  intergrid  operators  as  in  Case  1,  the  (1,  l)-component  of 
the  matrix  symbol  of  the  coarse  grid  corrector  CG  (see  (7))  for  low  frequen¬ 
cies,  denoted  as 

9^  —  a  =  (01,0:2)  w  0, 
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can  be  approximated  by 


CG{e,z)n  «  1- 


•2  +  Lff{2a) 


7/3(q:i  +  +  i{{asin{2ai)  +  bsin{2a2))  12  —  {asinai  +  bsina2)) 

h{z  +  Lh{2o.)) 

^  7^(of  +  -  i{aai  +  ba2)/2 

""  hiz  +  LH{2a))  ■ 

Taking  — ^  to  be  the  imaginary  part  of  Lf{{2a)^ 


CG{e,z)n  ~ 


7/?(ai  +  0:2)^  —  i{Q,a\  + 
8^(af  +  a|)2 


which  diverges  to  00  as  a  0  when  aa^  +  ba^  7^  0.  However,  for  the  related 
static  two  grid  iteration,  2  =  0,  one  has 

_ 7^(a^  +  Oilf  -  i{aal  +  bal)/2 _ 

^  ~  S/3{al  +  a^y  +  i{{aa\  +  602)  -  2(aaf  +  )/3)’ 

As  a  0,  it  converges  to  0  except  in  the  direction  of  aai  +  602  =  0.  In 
this  direction,  CG{6^z)ii  converges  to  either  7/8  or  6/8  depending  upon 
whether  aa\  +  ba\  vanishes. 

Since  the  effectiveness  of  a  coarse  grid  corrector  is  strongly  influenced  by 
its  action  on  the  error  that  has  very  low  frequencies,  the  above  observation 
indicates  that,  instead  of  reducing  the  error  in  smooth  modes,  the  coarse 
grid  corrector  in  the  two  grid  WR  iteration  could  magnify  the  error,  causing 
a  divergence  of  the  iterative  process.  Meanwhile,  its  static  analogue  works 
weU  [5]. 


4  Numerical  Results  and  Conclusions 

The  numerical  experiments  were  carried  out  on  two  case  studies  discussed 
in  the  previous  section.  The  problems  were  solved  by  the  two  grid  Jacobi 
WR  method  over  the  unit  square  along  0  <  t  <  tf,  tf  =  1  (it  has  been 
found  that  the  dependency  of  the  measured  convergence  rate  on  the  length 
of  time  interval  can  be  ignored  [16]).  The  space  derivatives  were  discretized 
as  described  in  Section  3  with  the  uniform  fine  grid  size  h  =  1/64.  The  initial 
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guess  was  randomly  generated  to  excite  all  possible  Fourier  modes.  The 
trapezoidal  rule  was  used  with  step  size  ,01  for  the  time  integration.  In  both 
cases,  the  asymptotic  ratios  of  the  defect  (see  (3)) 

i  =  (16) 

maxjg[o,i/]||4'  ^^(011 

were  collected  and  their  mean  value,  denoted  by  f,  was  used  as  an  approxi¬ 
mation  to  the  spectral  radius  p{V). 

Tables  3  and  4  list  the  range  of  the  ratios  for  Case  1  and  2  respec¬ 
tively.  Table  5  lists  the  defects  coUected  from  the  experiments  in  Case  2 
with  a  =  6  =  1  and  /?  =  .125,  using  and  without  using  the  multigrid  tech¬ 
nique.  The  results  show  that,  in  Case  1,  the  multigrid  technique  accelerates 
the  convergence  of  the  WR  iteration  in  a  similar  way  as  it  does  on  static 
iteration  for  the  corresponding  steady  state  problem.  However,  in  Case  2, 
the  multigrid  WR  process  quickly  diverges.  The  results  in  Table  5  clearly 
indicate  that  it  was  the  coarse  grid  correction  that  caused  the  divergence, 
as  shown  by  the  previous  analysis. 

All  the  numerical  data  confirm  the  results  obtained  by  the  Fourier- 
Laplace  analysis.  This  demonstrates  that,  for  a  converging  iterative  pro¬ 
cess,  the  analysis  provided  practical  estimates  to  the  actual  convergence 
rates;  while  for  a  diverging  iterative  process,  it  was  able  to  locate  the  source 
of  divergence,  which  can  be  very  helpful  in  searching  new  directions  and 
developing  better  methods  for  accelerating  the  WR  iteration. 
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Table  2:  Spectral  Radius  of  Two  Grid  Red-Black  Gauss-Seidel  WR  (6  =  1) 


a 

P  =  A 

13  =  .5 

j8  =  .75 

/?  =  !. 

1.0 

1.0  (1.0) 

.5189  (.5189) 

.4999  (.4999) 

0.5 

.8789  (.8789) 

.5960  (..5960) 

.4938  (.4882) 

.4914  (.4886) 

0.0 

.6574  (.6339) 

.5373  (.5133) 

.4967  (.4257) 

.4947  (.4119) 

Table  3:  Case  1  (w  =  2/3,  a  =  6  =  1) 


■■ 

f 

p(V) 

■iMil 

.6583 

.7086 

0.75 

.5349  -  - .7090 

.6354 

Table  4:  Case  2  (w  =  1/2,  a  =  6  =  1) 


/? 

r(2) 

rW 

”7^ 

r(^) 

,.(6) 

“777“ 

-757- 

7.(10; 

0.125 

0.9 

0.5 

1.3 

4.1 

3.6 

4.4 

3.7 

3.0 

2.8 

2.3 

0.5 

0.7 

0.4 

0.7 

1.3 

2,8 

3.0 

2.8 

2.4 

2.0 

2.2 

Table  5:  Case  2:  Measured  Defect 


Method 

2  =  3  i  =  4  ?  =  5  z  =  6  2  =  7 

Two  Grid  WR 

70.2  286.5  1017.7  4521.0  16512.6 

WR 

62.6  47.9  41.3  38.7  37.7 
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